Planetesimal and Protoplanet Dynamics in a Turbulent Protoplanetary Disk: 

Ideal Unstratified Disks 

Chao-Chin Yang 

Department of Astronomy, University of Illinois, Urbana, IL 61801 
Department of Astrophysics, American Museum of Natural History, New York, NY 10024 

cyangOamnh . org 
Mordecai-Mark Mac Low 
Department of Astrophysics, American Museum of Natural History, New York, NY 10024 

mordecai@amnh . org 
and 

Kristen Menou 

Department of Astronomy, Columbia University, New York, NY 10027 
kr istenOastro . Columbia . edu 

ABSTRACT 

The dynamics of planetesimals and planetary cores may be strongly influenced by 
density perturbations driven by magneto-rotational turbulence in their natal proto- 
planetary gas disks. Using the local shearing box approximation, we perform numerical 
simulations of planetesimals moving as massless particles in a turbulent, magnetized, 
unstratified gas disk. Our fiducial disk model shows turbulent accretion characterized 
by a Shakura-Sunyaev viscosity parameter of a ~ 10~ 2 , with root-mean-square density 
perturbations of ~10%. We measure the statistical evolution of particle orbital prop- 
erties in our simulations including mean radius, eccentricity, and velocity dispersion. 
We confirm random walk growth in time of all three properties, the first time that 
this has been done with direct orbital integration in a local model. We find that the 
growth rate increases with the box size used at least up to boxes of eight scale heights 
in horizontal size. However, even our largest boxes show velocity dispersions sufficiently 
low that collisional destruction of planetesimals should be unimportant in the inner 
disk throughout its lifetime. Our direct integrations agree with earlier torque measure- 
ments showing that type I migration dominates over diffusive migration by stochastic 
torques for most objects in the planetary core and terrestrial planet mass range. Dif- 
fusive migration remains important for objects in the mass range of kilometer-sized 
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planetesimals. Discrepancies in the derived magnitude of turbulence between local and 
global simulations of magneto-rotationally unstable disks remains an open issue, with 
important consequences for planet formation scenarios. 

Subject headings: Accretion, accretion disks — methods: numerical — MHD — plane- 
tary systems: formation — planetary systems: protoplanetary disks — turbulence 



1. INTRODUCTION 



The core accretion scenario is a key ingredient of ou r current theory for planet formation 
(see reviews bv lLissauerlll993l ; iPapaloizou Terquemll2006l . and references therein). In this model, 
solid bodies ranging from micron-sized dust grains to kilometer-sized planetesimals are the building 
blocks from which planetary cores assemble and ultimately form giant planets orbiting stars. The 
dynamics of planetesimals and planetary cores may be strongly influenced by their natal protoplan- 
etary gas disks. In particular, type I migration due to angular momentum exchange with the gas 



disk c ould drive a planetary core into the host star in a short time scale (e.g.. lGoldreich Sz Tremaine 



1980 



Ward 



1997; 



Tanaka et al 



2002 



Hartmann et al 



Menou &: Goodman! 120041) compa r ed to a typical di sk life 



1998; ISicilia-Aguilar et al. 



2006 



Hillenbrand! I2008T ). This 



time of <10 Myr (e.g., 

remains one of the major obstacles to planet formation in the core accretion scenario unless re- 
cent models of the countervailing effects of type III migration in non-adiabatic disks ar e confirmed 
(|Paardekooper fc Mellemal boOfil : klev fc Cridall200» baardekooper fc Papaloizoul boOfll ) . 



Turbulent transport in young protoplanetary disks seems required for mass accretion to occur 
at observed rates. The most likely caus e of this turbulence is tha t disks are unstable to linear pertur- 

Although this magneto-rotational 



bations if they are weakly magnetized (jBalbus &: Hawlev 



1991 



instability (MRI) is the most promising mechanism to drive turbulence in protoplanetary disks, 
the numerical convergence of simulations with increasing resolution continues to be an open issue. 
It has been shown that MRI-driven turbulence in ideal, magnetized, unstratified disks without 
mean magnetic fields decreases with increasing resolution and might be negligible when resolu - 
tion is high and numerical dissipation is small (jFromang &: Papaloizoul 120071 : iPessah et al.1 120071 ) . 
To maintain a nonzero, convergent le yel of turbulent accretion, it ha s been argued that simu- 
lated disks might require stratification (IDavis et al.ll2009l : IShi et al.ll2009l ). explici t dissipation (e.g. 



Simon et al 



Fromang et al.l 120071 : iLesur k, Longarettil 120071 ) . or nonzero magnetic flux (e.g., iGuan et al.ll2009l : 



2009h . Furthermore, the conductivity near the mid-plane of a protoplanetary disk may 



be too low for the MRI to o perate, at least over some range of radii, resulting in a quiescent region 



known as a dead zone (e.g. . iGammie 



1996 



Sano et al.l 120001 ; iFromang et al.l [2002J; ISemenov et al 



2004 : lllgner &; Nelson! 120061 : 1 Turner et al.l 120071 ) . Disks containing dead zones will, nevertheless, b e 
stirred by turbulent regions near the disk surface ([Fleming Stonell2003l : lOishi &: Mac Lowl l2009). 



The turbulent nature of protoplanetary disks might lead to scenarios for the formation and 
migration of planetesimals and planetary cores that are rather different from those suggested for the 
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better-studied case of a laminar disk. The turbulence driven by the MRI causes density enhance- 
ments significant enough to exert gravitational t orques that turn the orbital motion of low-m ass 
objects into a random walk (jLaughlin et alj|2004l . hereafter LSA04; iNelson fc Papaloizoul 12004 ) . It 
has been shown that radial exc ursions and e ccentricities of planetesimals or pro toplanets could be 
excited through such a process (|Nelsonll2005l . hereafter N05; lQgihara et alj|2007l . hereafter OIM07). 
Since some protoplanets could diffuse their way radially outward in the process, type I migration 
might be effectively d elayed and some ob jects might more easily survive past the gas disk depletion 



(Johnson et al 



2006, hereafter JGM06; lAdams fc Blochl 12003 ). Even if the disk contains a dead 



zone, lOishi et al.l (|2007l . hereafter OMM07) have shown that low-mass objects within th e dead 
zone still experience some of the turbulent torques generated by the active layers. Ilda et al.l (|2008l . 
hereafter IGM08) suggest, on the other hand, that the velocity dispersion of planetesimals excited 
by hydromagnetic turbulence might be so strong that kilometer-sized objects suffer from collisional 
destruction. 

The survivability of planetesimals or protoplanets under type I migration or collisional destruc- 
tion sensitively depends on their orbital dynamics in a turbulent gas disk. Previous direct orbital 
integrations of planetesimals or protoplanets embedded in a turbulent gas disk were conducted 
in global disk models (LSA04; N05; OIM07; IGM08). In contrast to global disk models, it can 
be advantageous to employ the local shearing box approximation (e.g., iGoldreich &: Lvnden-Bell 



1965 



Hawlev et al 



1995 



Brandenburg et al.lll995l ) because of it s high resolving powe r on t urbu 



lence structures and the possibility of integrating for long times. INelson fc Papaloizoul (|2004l ) and 
OMM07 first measured the stochastic torques generated by hydromagnetic turbulence at the center 
of a local shearing box. In this paper, we pursue this topic further by using direct orbital integration 
of planetesimals moving as massless particles under the gravitational influence of MRI-driven tur- 
bulence in a local shearing box. We focus our attention on unstratified disks in the context of ideal 
magnetohydrodynamics (MHD). To maintain a nonzero, numerically convergent level of stochastic 
perturbations driven by the MRI, we impose a constant vertical magnetic flux. We describe our 
numerical models in $2] and present the simulations in along with statistical analyses of the disk 
properties and the planetesimal orbits. In $31 we use our results to revisit the issue of survivability 
of planetesimals and planetary cores, before reaching our conclusions in $5) 



NUMERICAL MODELING 



We use the parallelized, cache-efficient, Pencil Code described by lBrandenburg k, Doblerl (J2002J) 
It solves the non-ideal MHD equations by sixth-order finite differences in space and third-order 
Runge-Kutta steps in time. The induction equation is solved using the magnetic vector potential 
A so that zero divergence of magnetic field B is guaranteed at all ti me. To save mem ory usage, the 
Runge-Kutta time integration is performed using the 2A^- method (| Williamsonl 1 1980 ) . The scheme 
is not written in conservative form. Instead, conserved quantities like total mass are monitored 
to evaluate the quality of the solution. In the following subsections, we describe the equations 
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assumed in our models as well as the numerical constructs. 



2.1. Magnetohydrodynamics 



We use the local shea ring box approximation (e.g., Goldreich Lynden-Bell 1965 : Brandenburg et al 



1995 



Hawley et al.lll995l ) to simulate a small Cartesian box carved out of a Keplerian disk at a 



large distance from the host star. The center of the box co-rotates with the disk at Keplerian an- 
gular speed Qk, the x-axis is directed radially, and the y-axis is directed azimuthally. The vertical 
component of gravity from the host star is ignored and thus the disk is unstratified. We impose 
a vertical, external magnetic field B cxt = B cxt z to maintain a finite magnetic flux. The MHD 
equations then become 



d t p - -Q K xd y p + V • (pu) 
3 

dtu £Ik x 9 v u + u ■ Vu 



d t A - -Q K xd y A 



(^Vp) , 

1. 



/d + V 

1 

P 



1 



+ - J x 

P 



B + B 



ext ) 



+fv + -V {v s pV • u) , 

p 

3 

-n K A y x + it x (B + B cxt ) + f R - rj s J, 



(1) 

(2) 
(3) 



in which p is gas density, u is gas velocity relative to the background shear flow, p is gas pressure, 
J = V x B / po is the electric current density, B = V x A, and po is permeability. The terms fo, 
f v , and f R , and those containing scalar variables v a and r] s are numerical dissipation terms needed 
to stabilize the scheme, which are described below. They are needed to resolve shocks, and because 
the difference scheme formally has vanishing dissipation. We assume an isothermal equation of 
state, p = c 2 s p, where c s is the isothermal speed of sound. 

The mass diffusion term V • (za^V p) in the continuity equation (pQ) and the bulk viscosity term 
V (vspV ■ u) I p in the momentum equation ([2|) are i mplemented to broaden shocks. The artificial 



kinematic viscosity v s is of von Neumann type (c.f., lHaugen et al. 



2004 



Lvra et al.ll2008h : 



0, 



h 2 V ■ u, if V • u < -c s /4h, 
otherwise, 



(4) 



where h is grid spacing. It is smoothed by taking a maximum over nearest neighbors and then 
convolved with a Gaussian kernel having a standard deviation of h. Note that the threshold for the 
velocity divergence is set for eliminating artificial diffusion where hydrodynamic shocks are unlikely 
to be present. 

We also include the Ohmic term —r\ s J in the induction equation fl3l ) to broaden st rong current 
sheets. The artificial resistivity r\ s assumes the same form used by iNitta et al.l (|200ll ) but with a 
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lower cutoff: 



HohvA, if v A > 8c s , 
0, otherwise, 



(5) 



where va = \B + 5 cxt | / yjpop is the Alven speed. Although simple, this form may perform better 
to resolve sha rp magnetic struct ures than a resistivity proportional to the magnitude of current 
density itself ([Fragile et al.ll2005l ). We apply a lower cutoff in equation ([5]) to only treat regions 
where fast magnetic reconnection may occur. 

In addition to applying artificial diffusions addressing shocks and current sheets, we also im- 
plement hyper-diffusio n fp, hyper- viscosity f v . and hyper-resistivity / p in t he respective MHD 



2004 



equations (P)-© (e.g.. lHaugen fc Brandenburg 

Sd = ^ V 6 p, fv = ^ ( V 6 « + S • V In p 



Johansen &: Klahr 



20051 ): 



and f R = i/ 3 V 6 A, 



(6) 



where is a constant, the tensor S is defined by Sij = djUi, and the sixth-order differential 
operator V 6 = + + d®. These terms are included in order to stabilize the high-order finite- 
difference scheme implemented by the Pencil Code. The corresponding diffusivity is proportional 
to k 6 , where k is the wavenumber of a signal in the simulation, so features at small scales dissipate 
at a much higher rate than those at larger scales. Therefore, by adjusting the coefficient V3, we 
can maintain numerical stability by damping oscillations near the grid scale while still preserving 
much of the inertial range of the modeled turbulence resolved in the simulation. In our models, we 
choose V3 such that the mesh Reynolds number 



Re 



mesh 



7T 



<1, 



(7) 



where u max is the absolute maximum of u over the computational domain. This criterion states 
that dissipation at the Nyquist frequency should be comparable to or stronger than gas advection. 
Notice that since all the adopted hyper-diffusive terms have the same coefficient 1^3, the effective 
magnetic Prandtl number PrM, c ff = Mo^eff / '"Heft in our simulations should be reasonably close to 
unity, where v e s and r) e g are the effective viscosity and resistivity in the simulations, respectively. 



2.2. Particle Dynamics 

In this work, we consider particles of zero mass to study the effect of hydromagnetic turbulence 
on the orbital properties of planetesimals. Although simple, this offers a good approximation for 
kilometer-sized planetesimals as they are large enough that gaseous drag force can be neglected, 
but small enough that type I migration does not dominate. For the approximation to be valid, 
the mass of a particle must lie between roughly 10 14 g and 10 26 g, corresponding to a size range 
of about 0.1-1000 km (OMM07). Furthermore, if the action of hydromagnetic turbulence on the 
particles is separable from effects due to other interactions, our measurements can be applied to all 
stages of planet formation. 
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In our models, therefore, we consider particles moving under only the gravitational influence of 
the host star and that of the protoplanetary gas disk. We ignore drag forces between the particles 
and the gas and the gravity of the particles. Because the particles exert no force on the gas 
or themselves, no migration torques act. Given deterministic Keplerian shear flow and epicycle 
motions, deviations in particle trajectories due to turbulent fluctuations in the gas can easily be 
isolated. 



Under these assumptions, the equations of motion for each particle become 



dx p 

~df 
du p 

~df 



Up - -Q K x p y, 



(8a) 
(8b) 



The vector x p is the position of the particle in the shearing box, while u p is the velocity of the 
particle relative to the background shear flow. The scala r variable $ is the gravitational potential 
of the gas, which is the solution of the Poisson equation (jGoldreich Lynden-Belllll965l ): 



d x -^Q K ta y ) +d 2 y + d 2 z 



(9) 



where G is the gravitational constant. 



Sheared, periodic boundary conditions (jHawlev et al.l Il995l ) are adopted for the solution of 
equation ([9]). The system is strictly periodic in the y- and z-directions while the x-direction requires 
special treatment. Since the yz-plane at any given x moves with the Keplerian shear flow, the lower 
x-boundary plane should be imaged by the upper x-boundary plane shifted by —3flKL x 5t/2 in the 
y-direction, where 5t is the time-step and L x is the x-dimension of the computational domain. 
In mathematical notation, p(x,y,z) = p(x + L x ,y — 3£lKL x 5t/2, z). Rather than interpolating 
p(x, y, z) in real space to obtain sheared periodicity, we use Fourier interpolation when solving 
equation Q in Fourier space (see Johansen et al. 2007 . Supplementary Information) 



The position x p and velocity u p of each particle is updated by solving the equations of mo- 
tion (|8|) simultaneously with the third-order Runge-Kutta steps for the MHD equations (d])-©. In 
addition to the Courant conditions set by the MHD equations, the time-step is limited by the abso- 
lute maximum of equation (|8ap such that no particles can cross more than half the zone size in one 
time-step. We compute the gradient of the potential V$ on the grid after solving equation ([9]) and 
then quadratically interpolate it to the position of each particle in the calculation of equation (|8bj) . 



lr This technique was originally suggested by Colin McNally at http: //imp. monaster . ca//,7Ecolinm/ism/rotf f t .html. 
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2.3. Code Units and Scaling Relations 

We define the length and the time units as the vertical scale height H and the orbital period 
P = 2h/Q,k, respectively, at the center of our local shearing box located at an arbitrary distance 
to the host star R. Since vertical hydrostatic equilibrium of isothermal gas requires that H = 
y2cg/Qjc, the speed of sound is fixed at c s = tt\^2. Note that this choice makes the system 
invariant with temperature. 

We adopt two different mass units such that po = (AirGP 2 )~~ 1 and po = (GP 2 )" 1 for a low- 
mass and a high-mass disk, respectively, where po is the uniform initial gas density. For these two 
disk models, the Toomre Q parameter for the gas is Q„ = c s VLk /ttGT, = 63 and 5.0, respectively, 
where S = ^ Po H is the column density! The gas dl sks in our models are gravitationally stable 
and thus we ignore gas self-gravity. For convenience, we define a dimensionless parameter 

£ = AttGpqP 2 = 4(2iT) 3 / 2 /Q g (10) 

as a measure of the strength of disk gravity. For our low-mass and high-mass disks, £ = 1 and 4n, 
respectively. 

In physical units, po is given by the following scaling relation: 

po = (1.2 x ID" 9 g cm" 3 ) £ (£) ~ 2 = (1.2 x 10 -9 g cm"*) t (A) ~* , (11) 

where M* is the mass of the host star. The corresponding column density is 

These scaling relations describe families of disk models to which our results apply. In particular, at 
1 AU around a solar-type star, the column density of our low-mass disk is roughly consistent with 



that of the classical minimum mass solar nebula (MMSN; lHavashilll98ll ). 

Finally, the units of magnetic field and vector potential are p}J 2 p\j 2 HP^ 1 and p]^ 2 p 1 ^ 2 H 2 P~ l , 
respectively. We arbitrarily set the permeability po = 1. The magnetic energy density associated 
with B = 1 is then l/(4-7r 2 ) of the initial pressure po = c 2 Po- 



2.4. Initial and Boundary Conditions 

The initial conditions for our models are the following. The gas density is uniform (p = po) 
while the magnetic vector potential is set to zero (A = 0). An external vertical magnetic field is 



2 Strictly speaking, this relation only holds for stratified disks. Comparison between stratified and unstratified disk 
models will be made in a subsequent study. 
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imposed, ranging from B ext = 0.01 to 0.64. The corresponding plasma (3 = 2p J Qpc 2 s / B 2 ranges from 
3.9 x 10 5 to 96. Gaussian noise in gas velocity of amplitude 10 3 is imposed to seed the MRI. 

As noted in §2.21 the boundary conditions for all dynamical variables are sheared periodic 
and we find values for the ghost zones using Fourier interpolation. Our fiducial model has a 
computational domain of 2 x 2 x 2H, but we also study domains with sizes up to 8 x 8 x 2H. The 
highest resolution we use is 64 grid points per disk scale height H. 

The coefficient 1/3 of the hyper-diffusive terms discussed in §2.11 needs to be fine-tuned such 
that the mesh Reynolds number Re mes h is as close to unity as possible during the course of the 
simulation (eq. [7]). We adopt an iterative approach to determine the optimal value of ^3, using 
Re mes h ~ 1- For the case of B cxt = 0.08 (/3 ex t = 6.2 x 10 3 ) with a 2 x 2 x 2H box at a resolution 
of 64 points per scale height, we choose = 2.9 x 10~ n for which u max ~ 6.0 ± 2.7 at saturation 
level, where the deviation of u mSuX is given by 3cr in its time variation. 

We uniformly distribute 32 3 particles in the entire computational domain. We do not allow 
them to move until time t = to after which the hydromagnetic turbulence has saturated and 
approached a statistically steady state. For the case of B cxt = 0.08 (f3 ext = 6.2 x 10 3 ), we choose 
to = 20P (see §3.1h . Then the particles are set to initially move relative to the background shear 
flow such that they have an initial eccentricity of eo and start at the apogee of their orbits, i.e., 

«p,o = -\HSl K (^_J y (13) 

(see Appendix lAl). We wrap a particle around when it moves beyond any of the six boundary 
planes. 

We remark that each model presented in the following section is just one realization of the 
stochastic nature of the turbulence, corresponding to one set of initial velocity perturbations of the 
gas. The similarity in particle orbital evolutions found across several models are due to closeness of 
the random number sequences used to generate the velocity perturbations and thus similar initial 
conditions for the gas. 



3. SIMULATION RESULTS 
3.1. Convergence of Turbulence Properties 

3.1.1. Grid Resolution 

We first present a study of the convergence with increasing numerical resolution of the prop- 
erties of the turbulence important to our work. We work on a 2 x 2 x 2H grid for this study to 
reach maximum resolution . Figure [1] plots density perturbation Ap/po, inverse plasma j3, and the 



Shakura fc Sunvaevl (|1973l ) a-parameter as a function of time t for disks with an external magnetic 
field of B ext = 0.08 ((3 ex t = 6.2 x 10 3 ) at resolutions up to 64 grid points per disk scale height 
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H, where Ap = p — p$. The density perturbation Ap/ pq shown is the root-mean-square value 
over the computational domain while the inverse plasma (3 shown is the volume-averaged value. 
The a-param eter is calculated fr om the combined effects due to the Reynolds and Maxwell shear 
stresses (e.g.. lBrandenburelll998l ): 



a 



y/2 {pu x u y - B x B y /p l0 ) 



P0C 2 s 



(14) 



where the bracket () denotes the volume average over the entire computational domain. As shown 
in Figure [TJ the MRI saturates and remains roughly steady after about t = 20P. After saturation, 
all three properties exhibit only small changes with increasing resolution, aside from a slight trend 
of increasing a. The properties of the saturated tur bulence appear to converge to a nonzero level, 
as opposed to disks without net magnetic flux (e.g., iFromang &: Papaloizoul 120071 ). 



3.1.2. Box Size 

We also examine how the same averaged quantities depend on the size of our local simulation 
box by studying three runs done with increasing horizontal size, up to 8 x 8 x 2H, at our medium 



resolution of 32 points per scale height. iJohansen et al.l (|2009l ) ran high-resolution (with ~137 points 



per scale hight), unstratified models without mean field and found roughly linear growth in the 
effective viscosity a with box size. In Figure [2] we show that the plasma (3 and the effective a- 
parameter appear almost independent of box size in our models, aside from a trend toward reduced 
temporal fluctuations when averaged over larger boxes. The root-mean-square density perturbation 
Ap/po shows a weak trend towards increasing at larger box size, though, with the average over the 
time period 20 < t/P < 120 increasing by 26% from the smallest to the largest box, a scale change 
of a factor of four. This probably occurs because of the inclusion of larger-scale instability modes 
in the larger boxes; it is not a particularly dramatic effect, though, because the smallest box size 
that we study already captures the fastest growing modes. Nevertheless, this small effect seems to 
strongly affect particle orbital properties, as we will discuss below. 



3.2. Vertical Net Flux Dependence 



We next turn to the effect of varying external vertical magnetic field. Figure [3] plots the same 
properties of the saturated turbulence as a function of the external vertical magnetic field, repre- 
sented by the inverse plasma f3, for our smallest box. They are volume averaged as described above, 
and then time-averaged over a period of at least 20P after saturation. Also included in the figure 
are the time variation of these properties, as indicated by the error bars. We confir m the general 



trend of increasing tu r bulence activity with increasing uniform vertical field (e.g., lHawlev et al 



19951 ; ISano et al.ll2004l ; IJohansen et al.l 120061 ) . Numerical convergence can be seen over the range 
of the field strengths we have explored. In addition to the turbulent transport often discussed in 




tip 

Fig. 1. — Density perturbation Ap/po, inverse plasma /3, and a-parameter as a function of time 
t (in units of orbital period P) for a 2 x 2 x 2H local shearing box under an external vertical 
magnetic field of -B ex t = 0.08 (/3 ex t = 6.2 x 10 3 ). All properties are volume-averaged over the whole 
computational domain with the root-mean-square value of Ap/po being given. Results are shown 
for resolutions up to 64 points per scale height H. 




tip 

Fig. 2. — Density perturbation Ap/po, inverse plasma j3, and a-parameter as a function of time t 
for three different box sizes at a resolution of 32 points per scale height H. An external vertical 
magnetic field of B ext = 0.08 (/3 cx t = 6.2 x 10 3 ) is imposed. Properties are volume-averaged over 
the whole computational domain and the root-mean-square value for Ap/pQ is shown. 
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the literature, we also report the dependency between density perturbation and external field in 
Figure O which may be more relevant to the orbital dynamics of particles moving in these disks. We 
emphasize that by varying the net vertical magnetic flux through a disk, a wide range of turbulent 
viscosity values can be obtained, as suggested by numerous previous works as well as Figure [3l 

To best represent typical protoplanetary accretion disks, we adopt a fiducial disk model with 
Bext = 0.08 (/? ox t = 6.2 x 10 3 ), which we run at a resolution of 64 points per scale height on a 2 x 2 x 
2H grid. As shown in Figures Q] and [31 this model gives a turbulent accretion of a ~ 10~ 2 , which 



is consisten t with current estimates for disks around typical T Tauri stars (e.g., lHartmann et al 



19981 . 120Q6I ) . Note that in our fiducial model, the root-mean-square density perturbation of the gas 



is on the order of ~10%. 



3.3. Motion of a Single Particle 

Using our fiducial model in the context of a protoplanetary disk, we now study the orbital 
dynamics of zero-mass particles moving in this turbulent environment. As demonstrated in Figured! 
a particle moving under gravity of the turbulent gas undergoes epicycle motion horizontally as well 
as continuous change in its mean radius. We define the radial drift of each particle as Ax = x — xq, 
where x is the mean radial position over one orbital period (as exemplified by the red line in Fig. 
and xq is the initial radial position. Given that x <C R, the eccentricity of each particle can be 
approximated by e « (x max — x m i n ) /2R, where x max and x m \ n are the maximum and the minimum 
radial positions in one epicycle, respectively. With these two quantities, we can measure the orbital 
migration and eccentricity change of planetesimals induced by hydromagnetic turbulence. 

Figure [5] shows the change of radial drift and eccentricity with time for four randomly selected 
particles. It is evidently a stochastic process and the final outcome can be quite different with 
slightly different initial conditions. Nevertheless, statistical methods can be employed to quan- 
tify the process. We discuss the statistical evolution of these orbital properties in the following 
subsections. 



3.4. Radial Drift 



Histograms of the distribution of radial drifts at three different times in the low-mass disk 
version of our fiducial model are plotted in Figure [H The distribution of particles in radial drift 
resembles a normal distribution with its center located at approximately zero. This is not surprising 
since there is no preferred direction locally for the turbulence to generate a net torque. More inter- 
estingly, the width of the distribution increases with time. Although the hydromagnetic turbulence 
has no net effect on the orbital radius of the particles, it becomes more and more likely for any 
single particle to drift away from its original orbit as time increa ses. This could help a subset of 
particles to survive type I migration, as suggested by JGM06 and lAdams fc Blochl ((2009). 




Fig. 3. — Density perturbation Ap/po, inverse plasma /3, and a-parameter as a function of external 
magnetic field in terms of inverse plasma (3 for a 2 x 2 x 2H box. All properties are volume-averaged 
over the whole computational domain as well as time-averaged over an interval of at least 20 orbital 
periods after saturation of the MRI. Results are shown for resolutions up to 64 points per scale 
height H. The error bars denote la in time variation around the volume-averaged properties. 
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Fig. 4. — Radial motion of one representative particle with initial eccentricity eo = moving in the 
low-mass disk version of our fiducial model. The black line shows its radial displacement x — xq 
from initial position x$ as a function of time t. The red line shows the corresponding radial drift, 
defined as the running average over one epicycle, i.e., one orbital period P. 



a: 
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Fig. 5. — Time evolution of radial drift Ax {top panel) and eccentricity e {bottom panel) for four 
randomly selected particles with initial eccentricity eo = moving in the low-mass disk version of 
our fiducial model. The eccentricity is in terms of H/R, the ratio of one disk scale height to the 
distance of the shearing box to the host star. 
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Ax/H 



Fig. 6. — Distribution of particle radial drifts Ax at three different times in our fiducial model. 
These particles have initial eccentricity bq = and evolve in a low-mass disk. The dotted lines are 
best fits in the form of normal distributions. 
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Figure [7] shows the standard deviation of radial drift a(Ax) as a function of elapsed time 
At = t — to for particles with different initial eccentricity eo moving in different strengths of 
disk gravity £ in our fiducial model. The standard deviation steadily increases with time, with 
little difference between particles of different initial eccentricities. Power-law fitting results in 
time indices of about 0.52-0.58, just slightly larger than 1/2. This confirms the proposition that 
gravitati onal influence of the hvdr omagnetic turbulence makes particles undergo random walks 



diffusion process (JGM06; OIM07; 



(LSA04: iNelson &; Papaloizoul 12004; N05) and the resulting: orbital evolution can b e described as a 



Adams Sz Bloch 



2009 



Rein fc Papaloizoul 12009). 



We measure the dependence of radial drift on disk gravity, as quantified by the dimensionless 
parameter £, by comparing the results from our low-mass and high-mass disk models (see §2.3|) . As 
demonstrated in Figure radial drifts in these models roughly coincide after being scaled with £, 
indicating a dependence close to linear. By assuming that a (Ax) scales with £Ai 1//2 , our best fit 
to the results shown in Figure [7] is 



a(Ax) = (3.8 ±0.4) x 10 -4 £ff 



At\ 



1/2 



(15) 



The size of the computational domain of a local shearing box does have a rather substantial 
effect on the magnitude of the random walk, and thus the diffusion derived from this model, as 
shown in Figure Increasing the horizontal size of the box by a factor of four results in almost an 
order of magnitude increase in the standard deviation of the radial drift at each time. In contrast 
to gas properties discussed in §3.1| we have not seen convergence of the magnitude of the orbital 
random walk with box size in our study. 

Our largest box shows an amplitude of the random walk roughly a factor of three smaller than 
seen in the global model described by N05. We make this estimate by examining his Figure 4. 
This shows the semi-major axis for six zero-mass particles at different radii R over a period of time 
of roughly 100 orbits in the relevant region of 2-3 times the inner edge of his grid, at R = R±, 
in a disk with H ~ 0.1R. The standard deviation of Ax/£H of the particles at At — 90P is 
roughly 0.06, where P is the orbital period at the respective initial radii of the particles and the 
relevant disk-gravity parameter £ ~ 0.9 (R/ R\) 2 (Richard P. Nelson 2009, private communication). 
By comparison, the result for our 8 x 8 x 2H box shown in Figure [5] gives a(Ax)/H ~ 0.02 at 
At ~ 90P. 

On the other hand, OIM07 found results roughly consistent with our fiducial disk model. 
These authors used orbital integrations of particles influenced by torques given by a heuristic, 
stochastic formula for hydromagnetic turbulence, and the formula was suggested by LSA04 based 
on zero net flux, global disk, MHD simulations. The disk models OIM07 studied were about 10— 
100 times less massive than the MMSN, but they reported their scaling with varying disk mass. 
By extrapolating their results to values appropriate for our low-mass disk model with £ = 1, and 
considering their fiducial magnitudes for the stochastic torques, we find that our measured spread 
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Fig. 7. — Standard deviation of radial drift a(Ax) as a function of elapsed time At in our fiducial 
model. The dotted lines are obtained from low-mass disks (£ = 1) while the dashed lines are from 
high-mass disks (£ = 4tt). Particles with initial eccentricities eo = 0, 0.1(H/ R), and 0.2(H/R) are 
denoted by red, green, and blue lines, respectively. The black solid line is the best fit to all six 
curves. 
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Fig. 8. — Standard deviation of radial drift a(Ax) as a function of elapsed time At for three 
different box sizes at a resolution of 32 points per scale height H (solid lines), where only the low- 
mass disk model (£ = 1) and particles with zero initial eccentricity are considered. For comparison, 
the straight dotted line gives the best fit to the high-resolution model shown in Figure [7] (eq. |15|). 



-20- 



of radial drift (eq. [15]) is roughly consistent with theirs at 1 AU. As discussed by these authors, 
large uncertainty in thier results might be involved, mostly due to the uncertainty in the magnitude 
of the stochastic torques and their neglecting the power given by LSA04 in the m = 1 mode (where 
the integer m represents the Fourier decomposition of density structure in azimuthal angle rather 
than the spiral mode used in density wave theories). 

3.5. Eccentricity 

Figure shows the histograms of eccentricity at three different times for particles with initial 
eccentricity eo = moving in the low-mass version of our fiducial disk model|f| The distributions 
appear close to a Rayleigh distribution: 




where a is a constant proportional to both the peak and the width of the distribution. The standard 
deviation of this distribution is given by (j\J (4 — 7r)/2. As shown in Figure [9^,, both the eccentricity 
at the peak value and the width of the distribution increase with time, so it at first appears that 
eccentricity is excited by hydromagnetic turbulence. Such eccentricity growth was also reported by 
N05 and OIM07. 

However, we find a different distribution for particles with nonzero initial eccentricity in the 
same model. Figure Eb shows histograms of eccentricity deviation Ae = e — eo at three different 
times for particles with initial eccentricity eo = 0.1(H/R). They are similar to a normal distribution 
with a constant mean, implying that the average eccentricity of these particles remains constant 
at the initial eccentricity (see also OIM07). The width of the distribution does increase with time. 
Therefore, hydromagnetic turbulence does not just excite the eccentricity of particles; it can also act 
to damp the existing eccentricity of some particles. The distribution found for particles with eo = 
is just a special case of this general behavior: since eccentricity is a positive definite quantity and the 
initial eccentricity is zero, it is no surprise that a normal distribution for the eccentricity deviation, 
/(Ae) = exp (-Ae 2 /2cr 2 ) /yfifta, manifests as a Rayleigh distribution for the eccentricity itself 
(eq. PE|). 

Figure [10] shows the standard deviation of eccentricity deviation cr(Ae) as a function of elapsed 
time At measured for our fiducial model. Note that to be consistent with a normal distribution as 
discussed above, we have multiplied the standard deviation measured for particles with eo = by a 
factor of -y/2/(4 — ir). As in the case of radial drift discussed in §3.41 little difference exists between 



3 We remark that the mean eccentricity of particles with initial eccentricity eo = is about 0.012(H/ R) and 
0.16(H / R) at t = 500P for the low-mass and the high-mass disk models, respectively. These eccentricities correspond 
to epicycle motions covering about 1.5 and 20 grid zones in the radial direction, enough to resolve the gravitational 
forces experienced by the particles. 





Ae/{H/R) 

Fig. 9. — Distribution of particles in (a) eccentricity e for particles with initial eccentricity eo = 
and (b) eccentricity deviation Ae = e — eo for particles with eo = 0.1(H/R) at three different times 
in our fiducial model. These particles move in a low-mass disk. The dotted lines are best fits in 
the form of (a) Rayleigh and (b) normal distributions. 
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Fig. 10. — Standard deviation of eccentricity deviation cr(Ae) as a function of elapsed time At in 
our fiducial model. The line styles and colors are the same as in Figure [71 The standard deviation 
measured for particles with initial eccentricity eo = has been multiplied by a factor of ^2/(4 — it) 
to be consistent with a normal distribution (see §3.5|) . 

particles with different initial eccentricities and the results scale linearly with the dimensionless 
parameter £, which measures the strength of the disk gravity. Chi-square fitting with a ^Ai 1 / 2 
dependence to the growth in eccentricity deviation in all six models leads to (Fig. fl"0|) 

a(Ae) = (4.1 ± 0.6) x 10" 4 £ (|) (^j ^ . (17) 

Note that at 1 AU in an MMSN disk, £ — 1, H/R ~ 0.1, and P = 1 yr, and thus the increase in 
eccentricity deviation due to hydromagnetic turbulence only amounts to about 0.04 in 1 Myr. 

As with the case of radial drift discussed in Q3A\ our measured spread of eccentricity in our 
fiducial model is in approximate agreement with OIM07. On the other hand, N05 reported a typical 
eccentricity growth of e ~ 0.03 for At ~ 100P. With £ ~ 6.6 at radius R = 2.7R\ in the model 
studied by N05 (Richard P. Nelson 2009, private communication) and H/R ~ 0.1, the eccentricity 
deviation given by equation (fT7]) is about one order of magnitude less than what was reported by 
N05. 

One possible explanation for this inconsistency could be the presence of large-scale structures, 
particularly m = 1 modes, in global, but not in local models. OIM07 reported that the inclusion of 





, n-4 I l l l l l l l l I l l l l l l l l 

IU 1 10 100 

At/ P 

Fig. 11. — Standard deviation of eccentricity deviation cr(Ae) as a function of elapsed time At 
for three different box sizes at a resolution of 32 points per scale height H (solid lines), where 
only the low-mass disk model (£ = 1) and particles with zero initial eccentricity are considered. 
For comparison, the straight dotted line gives the best fit to the high-resolution model shown in 
Figure Q3S (eq. [T7j). 



an m = 1 mode in their torque formula induces a ten times greater impact of hydromagnetic tur- 
bulence on particle orbits, using the calibration with global MHD simulations provided by LSA04. 
Inspection of LSA04, as well as N05, shows that the density structures in their global models often 
extend more than tt/2 in azimuthal angle, leading to the m = 1 mode having the largest ampli- 
tude. Such large-scale modes indeed can be excited in self-gravitating disks, and they could also 
be excited by local turbulence. This idea is offered some support by the modest growth in density 
pe rturbations in larger b oxes that we find (Fig. [2]), and the growth in a reported in larger boxes 
by Ijohansen et al.l (|2009l ) . 

In Figure [11] we examine the growth in eccentricity deviation as a function of box size in our 
local models. Increasing the horizontal box size by a factor of four indeed increases the eccentricity 
deviation by a factor of four, though this still does not account for the order of magnitude higher 
value found by N05 in his global model. If m = 1 modes indeed can be similarly excited by well 
resolved turbulence, though, this could offer an explanation for the discrepancy. 



However, using local models with mean azimuthal fields, 



Guan et al 



(|2009h reported that 
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the coherent structures induced by the MRI are localized, with correlation lengths of about 0.05-ff , 
0.32H, and 0.05H in radial, azimuthal, and vertical directions, respectively. Therefore, it is possible 
that current global models might not have enough resolution to model such fine structures, which 
were then spuriously connected into extended structures resembling m = 1 modes. To determine 
the physical reality of large-scale structure in turbulent disks, global models capable of resolving 
localized structures will be needed. 



Another important quantity in the study of planetesimal dynamics is the velocity dispersion 
of the particles. We calculate the radial and the azimuthal components of velocity dispersion by 
taking the standard deviations of radial and azimuthal velocities for all particles, i.e., cr(u x ) and 
a(u y ), respectively. Since all the dynamical equations are linearized in terms of x/R in the local 
shearing box approximation, the background shear flow is uniform irrespective of position, so the 
velocity dispersion for all particles in the computational domain is a well-defined local quantity. 
Figure [T2l plots cr(u x ) and cr(u y ) measured for particles with zero initial eccentricity in the low-mass 
disk version of our fiducial model as a function of time t. The velocity dispersion monotonically 
increases with time, and thus hydromagnetic turbulence tends to steadily heat up a planetesimal 
disk. Note that cr(u y ) ~ a(u x )/2, a sanity check that our results are consistent with a swarm of 
non-interacting particles moving epicyclically in a Keplerian disk. 

Figure [13] compiles the radial velocity dispersions <r{u x ) as a function of elapsed time At 
measured in our fiducial model for particles with different initial eccentricity eo moving in disks of 
different disk gravity. As is the case with the radial drift and eccentricity deviation discussed in 
§3.41 and §3.51 the results depend linearly on £, the strength of disk gravity, though there seems 
to be a slightly enhanced effect for particles with eo > 0. Given the uncertainty involved in the 
numerical simulations, we assume it is a secondary effect as a first approximation. The best fit to 
all six models is then given by 



The corresponding timescale for turbulent excitation of velocity dispersion can be estimated by 
r = fj/(da/dAt), and we find from equation (|18|) 



Increasing the box size increases the magnitude of this effect, as shown in Figure [T4"l possibly even 
non-linearly. 

As noted in £J33J hydromagnetic turbulence can act either to excite the eccentricities of plan- 
etesimals or to circularize their eccentric orbits. However, a perfectly cold disk of planetesimals 



3.6. Horizontal Velocity Dispersion 




(18) 




(19) 
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Fig. 12. — Radial and azimuthal components of velocity dispersion a(u x ) and cr(u y ) as a function 
of time t for particles with initial eccentricity eo = moving in the low-mass disk version of our 
fiducial model. They are presented as running averages over one orbital period P. 
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Fig. 13. — Radial component of velocity dispersion a{u x ) as a function of elapsed time At in our 
fiducial model, computed as running averages over one orbital period P. The dotted lines are 
obtained from low-mass disks (£ = 1) while the dashed lines are from high- mass disks (£ = 47r). 
Particles with initial eccentricities eo = 0, 0.1(H/R), and 0.2(H/R) are denoted by red, green, and 
blue lines, respectively. The solid black line is the best fit to all six curves, assuming no explicit 
dependence on eo- 
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Fig. 14. — Radial component of velocity dispersion a(u x ) as a function of elapsed time At for three 
different box sizes at a resolution of 32 points per scale height H {solid lines), where only the low- 
mass disk model (£ = 1) and particles with zero initial eccentricity are considered. For comparison, 
the straight dotted line gives the best fit to the high-resolution model shown in Figure [13] (eq. [18j). 
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with nonzero initial eccentricity but vanishing velocity dispersion will be monotonically heated 
by the turbulence, as shown above. The velocity dispersion increases with time while the mean 
eccentricity may remain unchanged unless a significant fraction of particles are circularized. There- 
fore, the eccentricity may not necessarily be proportional to the velocity dispersion for a swarm of 
planetesimals moving through hydromagnetic turbulence. 

The velocity dispersion of a planetesimal disk also grows with time due to mutual gravitational 
scattering. The c orresponding timescale for a swarm of identical particles of mass m p can be 
estimated by (e.g.. iPapaloizou fc Terquemll2006l ) 



a' 



T~GS 



(Ux) 



8y/7rG 2 m 2 n p In A p 




1 



(20) 



where n p is the number density of planetesimals and A p 



3a 2 ( 



)H p /4Gm p , in which H p 



\/2a{u z )/Q,K is the scale height of the planetesimal disk determined by the vertical velocity disper- 
sion a{u z ). To find tgs i n physical units, we assume for simplicity that most of the solid material 
in a protoplanetary disk is concentrated in planetesimals, and thus n p m p ~ epoc s /a(u z ) where e is 
the solid-to-gas ratio. We also assume a{u z ) ~ a(u x )/2, so that equation (f20|) becomes 



t~gs 



with 



Ae^c s Gm p In A p 



\/3 ln ( 2 + V3 
4 n \2-V3, 



iy p 



16irGm r , 



(21) 



(22) 



We further focus our discussion on a velocity scale of order v eBC , the escape velocity at the surface 
of a planetesimal: 



/32vr 



G s mi 



1/6 



(23) 



y 3 " " l pPp 

where p p is the material density of the planetesimal. This scale is of critical interest for planetary 
cores to accrete solid material; particles with relative velocities of order v eS c are more likely to 
coalesce into larger bodies than to be eroded into smaller pieces. By assuming a{u x ) ~ v esc , the 
timescales for heating a planetesimal disk by hydromagnetic turbulence and gravitational scattering 
become 

1/3 
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respectively. 

With the scales assumed above and at 1 AU in an MMSN disk, tqs ~ 2 x 10 4 yr while 
measurement from our high-resolution fiducial model described by equation (|19[) gives ~ 100 yr. 
For larger objects approaching the planetary mass regime, with m p = O.OOIMq, tqs ~ T ~ 
3 x 10 6 yr. Therefore, hydromagnetic turbulence probably dominates the heating of a disk of 
kilometer-sized planetesimals, while gravitational scattering may be more important for objects 
approaching Earth size. Models with larger boxes yield smaller values of tt, strengthening this 
conclusion. Note that tqs increases more steeply with a(u x ) and P but decreases less rapidly with £ 
than tt (eqs. [19] and |21j). and thus hydromagnetic turbulence gains dominance over gravitational 
scattering for larger velocity dispersion, larger distance to the host star, and more massive gas disks 
than assumed here. 



4. IMPLICATIONS FOR PLANET FORMATION 

In this section, we apply our results on the orbital evolution of zero-mass particles to two specific 
problems in planet formation. First, we estimate the strength of diffusive migration of protoplanets 
due to hydromagnetic turbulence and discuss their survivability under type I migration following 
the analytical framework established by JGM06. Secondly, we revisit the proposition of IGM08 
that planetesimals may suffer from collisional destruction as a result of hydromagnetic turbulence 
excitation of velocity dispersion among them. 



4.1. Diffusive Migration of Protoplanets 



Using a Fokker-Planck formalism, JGM06 derived an advection-diffusion equation to describe 
the evolution of the distribu tion of protoplanets un der the influence of both type I migration and 
hydromagnetic turbulence. lAdams k, Blochl (J2009J) further consolidated the analysis by studying 
more realistic disk density structure with both spatial and time dependence. These authors found 
that turbulence tends to reduce the lifetimes of most protoplanets while allowing some of them 
to linger long enough to survive rapid inward type I migration. The likelihood of producing a 
planetary system with a specific configuration sensitively d epends on the strength of the diffusive 
migration induced by the turbulence, however. JGM06 and lAdams & : Bloch (200 9j) calibrated th e 
turbulence strength with the global disk models computed by LSA04. INelson Papaloizoul (J200J), 
and N05. In contrast, using a local disk model, OMM07 found that the strength might be several 
orders of magnitude less than what was estimated by JGM06. In this section, we provide a new 
assessment based on the disk models studied in this work. 



As shown in §3.41 an initial delta function in the distribution of mean orbital radii of zero-mass 
particles is spread with time into a normal distribution of constant mean. A i 1 / 2 time dependence 
for the standard deviation of the distribution suggests that this process can be described by a 



- 30 - 



diffusion equation of the form 



where T> is the diffusion coefficient and / = f(t,x) is the distribution function: f(t, x)dx is the 
probability of finding a particle with a radial displacement in (x, x + dx) at time t. Let us write, 



/(M) = ^7^ exp 

1/2 



X 2 



(28) 



a(t) = (T1 V^pJ (t>0) ' (29) 

where o\ is a proportionality constant. Substituting equation (j28|) into equation ([27|) . we find that 
da/dt = IV I a and thus with a(t -> + ) = 0, 

<r(t) = 2P 1 / 2 t 1 / 2 . (30) 

By comparing equations ([29]) and ([30]) . the diffusion coefficient V is related to <ti and P by 

^ = 5" ( 31 ) 



The best-fit to our fiducial high resolution, small box model, given by equation (|15h . can be sub- 
stituted here to find 

£>(#) = 3.6 x 10~ 8 £ 2 (—J (32) 

by identifying t with At. In this derivation, we have assumed that the stochastic torque exerted 
by hydromagnetic turbulence is a local process such that the diffusion coefficient T> is sufficiently 
constant near x = 0. 

We are now in a position to estimate the diffusion coefficient D{J) of JGM06 in comparison to 

1 /2 

T>(R), where J = m p (GM+R) 1 is the orbital angular momentum of a protoplanet of mass m p or- 
biting a star of mass M* on a quasi-circular orbit at a radial distance R. Since D{ J) and T>(R) have 
dimensions of [J 2 /t] and [R 2 /t], respectively, they may be related by D(J) ~ V(R)(dJ/dR) 2 = 
(J /2R) 2 V(R). Using dimensional arguments, JGM06 defined a dimensionless parameter e to de- 
scribe the uncertainties associated with hydromagnetic turbulence: 

, n / o . . , o S 2 J 7 (2.1 x 10- 3 ) „ 2 fH\ 2 [J 2 \ . . 

D(J) = (2.1 x lO- 3 )^) 3 ^^ = '-^-^ee (-) (-) . (33) 



By comparing equations (|32j) and (l33]l . we find that 

e ~ 2.2 x 10~ 4 (34) 



for our fiducial model. Note that e is a constant independent of R given our scalings, in agreement 
with the assumption of constant e made by JGM06. 
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The value of e estimated in equation (|34p is about an order of magnitude smaller than what 
was reported by OMM07 for a stratified disk model with zero net flux on a 1 x 4 x AH grid 
at the same resolution as our fiducial model with 64 points per scale height. To identify the 
reason for this discrepancy, we further evaluate the magnitude and the correlation time of the 
torques exerted by the turbulent gas on the particles in our model. Firstly, the root-mean-square 
of the y-component of the gravitational force per unit mass exerted by the gas over all time and 
particles is a y , ims ~ (3.7 x 1(T 3 )£#P~ 2 = (4.2 x 10 -3 )(27rGE). We find negligible difference 
between torques calculated following the particles along their orbits and those calculated at the 
fixed center of the box. The magnitude we obtained is reasonably consistent with the value of 
Gj/.rms = (3.2 x 10 -3 )(27rGX) reported by OMM07 at the same numerical resolution. Secondly, we 
plot in Figure [15] the autocorrelation functions (ACF) of a y for several randomly selected particles 
as well as at the center of the box. The results are again similar to what was reported by OMM07, 
indicating a similar estimate for the correlation time, r c . The source of the discrepancy appears to 
be neither of these factors. 

We notice, though, that OMM07 could have overestimated the correlation time by equating it 
to the value of the second zero-crossing of the ACF. Strictly speaking, the correlation time should 
instead be computed by integrating an ensemble average of the ACF over all possible realizations, 
as motivated by the definitions of the diffusion coefficient and the correlation time in JGM06: 

D( J) = \J^ 5T(t - ~ J)6T(t + ~ J) dr (35) 

and 



t c = D(J)/5T 2 (t, J), (36) 

where 5T(t, J) is the fluctuating part of the torque and the overline denotes ensemble average. 
Since we distribute numerous particles uniformly over the entire computational domain, the ACF 
obtained for each particle can be considered as one realization, and the average of the ACFs over all 
particles may resemble the true ensemble average. This averaged ACF is shown by the solid black 
curve in Figure [T5l Note that the oscillation occurring at time lag longer than the second crossing 
for each particle is much reduced, indicating the noise nature of the autocorrelation at long time lag. 
However, the negative value of the ACF between the first two zero-crossings remains significant. 
This interval represents the anti-diffusion nature of the stochastic torques that we believe could be 
responsible for reducing the diffusion coefficient. Therefore, we suggest that a better approximation 
for the correlation time in accordance with equations ([35]) and ([36]) be 



f n °°ACF(r dr , , 

Jo V 1 — , 37 

2ACF(0) 



where ACF(r) is the ensemble-averaged ACF of the torque per unit mass a y as a function of time 
lag r. Using equation (|37p . we find that in our fiducial model r c ~ 0.020-P, about one order of 
magnitude less than the value of r c ~ 0.31P reported by OMM07. With the same approximation 
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Time Lag [P] 

Fig. 15. — Autocorrelation functions of the azimuthal (torque) component of the gravitational force 
per unit mass, a y , exerted by the turbulent gas in our fiducial model. The short-dashed black curve 
corresponds to a y calculated at the center of the box, while the various colored curves correspond 
to a y calculated following the orbits of several randomly selected particles. The average over all 
particles is given by the solid black line, which shows almost no power beyond the second zero- 
crossing. The vertical long-dashed black line indicates our estimated correlation time of the torques 
using equation ([37j) . Notice the negative part of the ensemble- averaged autocorrelation function 
beyond the first zero-crossing, which may be responsible for reducing the diffusion coefficient (see 
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D(J) ~ mpR 2 a,yT c used by OMM07, we find that e ~ 1.7 x 10 4 , in good agreement with the 
estimate of e ~ 2.2 x 10 -4 we derived from our direct measurement of particle radial drifts. 

Therefore, we have achieved consistent results using two independent approaches to estimating 
the diffusion coefficient, one by direct measurement as in equation (I32p and the other by the 
definition of correlation time given by JGM06 (eq. [36]). In principle, our direct measurement 
should be robust while using D{J) ~ t c 5T 2 can only be considered as an approximation. To 
use the latter approach, one should refer to equation (|37p to measure the correlation time of the 
stochastic torques, instead of conventional methods like using the zero-crossing of the ACF or the 
timescale of the peak of the temporal power spectrum, which probably give about one order of 
magnitude larger values. It will be enlightening for studies based on global models to conduct the 
same exercise described in this section, since this will likely explain part of the orders of magnitude 
discrepancy in diffusion coefficients derived from local versus global models. (Although the box size 
effects we have identified must also contribute to this discrepancy, they appear unable to increase 
the derived value of e by more than one order of magnitude.) 

According to Figures 6 and 7 of JGM06, the value of e we inferred from our simulations 
indicates that advective (type I) migration dominates over diffusive (stochastic) migration for the 
parameter space JGM06 have investigated. For an Earth-mass protoplanet at a radial distance 
up to 100 AU in the MMSN disk and in a viscous disk with a = 0.02, advection dominates when 
e < 10~ 2 and 10 _1 , respectively. For a protoplanet of mass as low as 0.01 M e at R = 10 AU in 
the same disk models, advection dominates when e < 10~ 3 and 10 -2 , respectively. The critical 
distance and mass for the transition between dominance of advection and diffusion given our small 
estimate of e for our fiducial model is outside of the parameter space explored by JGM06, and 
even our largest box has a value of e less than an order of magnitude larger. By inspection of 
Figure 7 in JGM06, though, the transition masses for an object at 10 AU in the MMSN disk and 
in the viscous disk probably lies at about 10 -3 and 10 -4 AT®, respectively. Therefore, our results 
suggest that hydromagnetic turbulence does not significantly affect the secular migration of Earth- 
sized protoplanets in regimes of current astrophysical interest. Nevertheless, torques exerted by 
turbulent density perturbations seem to be a dominant agent determining the orbital dynamics of 
kilometer-sized planetesimals. 



4.2. Collisional Destruction of Planetesimals 

IGM08 suggested that hydromagnetic turbulence may inhibit the growth of kilometer-sized 
planetesimals. They argued that the velocity dispersion excited by the turbulence could be so 
large that collisions between planetesimals exceed their material strength or self-gravity, leading to 
destruction. Their conclusion, however, relies on orbital integrations incorporating the heuristic, 
stochastic formulas for the time history of gravitational torques provided by LSA04, which in turn 
were calibrated using MHD simulations of a global disk model. As pointed out in §3.51 a possible 
inconsistency exists between global and local models. Since the latter shows a significantly lower 
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effect on orbital dynamics of planetesimals, it is worthwhile revisiting the planetesimal growth 
problem in light of the results highlighted in this work. 

For comparison purposes, we adopt the same scalable MMSN disk model as used by IGM08. 
The gas density and the speed of sound in the mid-plane are given by 

/ R \ ~ 11/4 

po = (1.8 x 1(T 9 g cm" 3 ) f g i^—j (38) 

and 

/ R \ ~ 1/4 

Cs = (l.lxl0 5 cms- 1 )!— J , (39) 

respectively, where a solar-type host star (M* = Mq) is assumed and f g is a scale factor. When 
f g = 1, the disk mass is about 1.4 times that of an MMSN disk. The corresponding ^-parameter 
and ratio of disk scale height to radial distance then become 



1/4 



e-^Hroj (40> 

and 

H ( R \ 1/4 , , 

— —0.051 (^j) • («> 

respectively. Substituting equations (|4"U|) and (|4"Tj) into equation (fTT)) . which is derived from our 
fiducial high-resolution model, we arrive at 



a(Ae)=3.1xlO- 5 / 9 (-^-) ( = ) . (42) 



(R_y l/A (^\ 1/2 



Comparing equation (|42p with equation (13) of IGM08 with the understanding that a(e) = a(Ae)y / (4 — ir) /2 
for particles with zero initial eccentricity (see §3.5[) , we find the value of the dimensionless param- 
eter 7 — a measure of the strength of hydromagnetic turbulence used by IGM08 — in our orbital 
integrations to be 7 ~ 2.0 x 10~ 4 . The rest of the analysis performed by IGM08 remains unchanged 
since the effects induced by hydromagnetic turbulence are all incorporated in the parameter 7. 

As noted in $\3.5\ enlarging the horizontal size of our shearing box by a factor of four at our 
medium resolution increases the amplitude of the eccentricity deviation by about a factor of four. 
Nevertheless, the largest box we have studied gives 7 ~ 6 x 10~ 4 , which remains somewhat smaller 
than the 7 ~ 10 _3 -10" 2 estimated by IGM08. 

The rather small values of 7 we obtained indicate that hydromagnetic turbulence might not 
pose as serious a threat to the growth of kilometer-sized planetesimals as suggested by IGM08. 
These authors compared the critical radii of planetesimals for accretive and erosive regimes due to 
different turbulence strengths 7 = 10 -2 , 10 -3 , and 10 -4 in their Figures 3, 4, and 5, respectively. 
Therefore, the values of 7 we measured point to a scenario in between what is predicted by Figures 4 
and 5 of IGM08. In this scenario, the erosive regime only appears in the outer regions of a 
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young protoplanetary disk and it disappears rapidly with decreasing disk mass. Kilometer-sized 
planetesimals may be able to evade collisional destruction in the inner regions of the disk. 

We reiterate that the value of 7 measured here pertains to a local region with turbulent 
stresses such that a ~ 10~ 2 and Ap/po ~ 10%. As discussed in §3.51 the validity of our results in a 
global context remains to be demonstrated by numerical experiments on a global disk model with 
a resolution capable of resolving the characteristic scales of coherent turbulence structures. 

5. CONCLUSIONS 

In this work, we have used local, shearing-box simulations to study the dynamics of massless 
planetesimals in a turbulent, isothermal, unstratified, gas disk driven by the MRI. With a uniform, 
vertical magnetic field but without explicit physical dissipation, the saturated turbulence is main- 
tained at a roughly constant level, showing convergence with increasing resolution. By adopting 
a suitable magnitude for the net magnetic flux, we produce a fiducial disk model with turbulent 
accretion at the level of a ~ 10 -2 and with root-mean-square density perturbations Ap/po ~ 10%. 
As discussed in §2.31 this model can be scaled to other physical systems of interest, as long as the 
assumption of negligible self-gravity of the gas remains valid. After the hydromagnetic turbulence 
in our fiducial model reaches saturation, we distribute numerous particles of zero mass and integrate 
their orbital motion under the gravitational influence of the turbulent gas. 

The stochastic nature of the orbital evolution of these particles is evident, so we characterize 
their orbital dynamics with statistical distributions, finding three major results. First, although the 
mean orbital radius does not change, particles slowly drift away from their original radii, so that 
the distribution of radii grows with time. Second, gravitational force from density perturbations 
produced in hydromagnetic turbulence can both excite and damp the eccentricities of particle orbits, 
with again no change in the mean value, but a growing width of distribution if the particles possess 
non-negligible initial eccentricities. Finally, the planetesimal disk is heated up by the turbulence, 
a process dominating over gravitational scattering between particles in most physical conditions 
relevant to protoplanetary gas disks and planetesimal sizes. A corollary of these results is that 
eccentricity does not serve as a good indicator of the velocity dispersion of the particles. 

The amplitude of orbital changes driven by the turbulence in our local models is significantly 
smaller than what was reported in recent global models (LSA04; N05; IGM08). Two possible 
explanations for this discrepancy suggest themselves: either insufficient resolution in the global 
models or the lack of convergence with box size in our local shearing box model, as discussed in 
§3.51 If our local results are valid, they indicate that although hydromagnetic turbulence can drive 
radial diffusion, eccentricity variations, and relative velocities of planetesimals and protoplanets, 
these effects may not be dominant in determining their evolution. In particular, it appears that 
type I migration dominates over turbulent radial drift for objects well above 10 -4 M®. In addition, 
hydromagnetic turbulence might not be exciting sufficient velocity dispersion to drive planetesimals 
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into an erosive regime that would inhibit their further growth. Before these results can be considered 
robust, however, it will be necessary to elucidate, and hopefully reconcile, the differences that have 
appeared between global and local models. 

We thank Jeffrey S. Oishi and Richard P. Nelson for the clarification of their works and useful 
discussions. We thank the anonymous referee for urging us to explore models with larger box 
sizes. The computations reported here were performed on the Columbia and Pleiades systems of 
the NASA Advanced Supercomputing (NAS) Division. This research is supported by the NASA 
Origins of Solar Systems Program under grant NNX07AI74G. 



A. VELOCITY OF A PARTICLE AT THE APOGEE OF ITS ORBIT 

In this section, we (re-) derive the velocity of a particle at the apogee of its elliptical orbit in 
the local shearing box approximation of a Keplerian disk. We repeat the equation of motion (I8bj) 
for a single particle without the gravity of the gas here: 

^ = 2Q K u y , (Al) 
^ = ~\n K u x , (A2) 



where we have dropped the subscript p for clarity. Eliminating u y in equations (|Aip and (|A2p leads 
to 

^ + n 2 K u x = 0. (A3) 

By assuming the particle is at the apogee at t = 0, the solution for u x is 

u x = -AsinQ^i, (A4) 

where A is the amplitude of the radial velocity. Since dx/dt = u x (eq. |8aj). the radial oscillation 
is then 

Ax = x — xq = — — cos {lxt, (A5) 
\l K 



where xo is the radial position of the center of the orbit. From equation (1A2|) . the corresponding 
azimuthal velocity relative to the background shear flow is 

u y = —-Acosttxt = — -^ft"Ax. (A6) 

Since the eccentricity e of the orbit is related with the amplitude of the radial oscillation by 
e « Ax(t = 0)/R, where R is the distance to the central object, 

««/(< = 0) = -\RV K e = —HSIk > ( A7 ) 

where H is the disk scale height. This is just the initial condition (|13p we set out to prove. Note 
that we have normalized e by the ratio H/R. 
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